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SUMMARY 

The various formulations of Maxwell’s equations are reviewed with emphasis on those 
formulations which most readily form analogies with Navier’s equations. Analogies involving 
scalar and vector potentials and electric and magnetic field components are presented. 
Formulations allowing for media with dielectric and conducting properties are emphasized. It 
is demonstrated that many problems in electromagnetism can be solved using the NASTRAN 
finite element code. 

Several fundamental problems involving time harmonic solutions of Maxwell’s equations with 
known analytic solutions are solved using NASTRAN to demonstrate convergence and mesh 
requirements. Mesh requirements are studied as a function of frequency, conductivity, and 
dielectric properties. 

Applications in both low frequency and high frequency are highlighted. The low frequency 
problems demonstrate the ability to solve problems involving media inhomogeneity and 
unbounded domains. The high frequency applications demonstrate the ability to handle 
problems with large boundary to wavelength ratios. 


INTRODUCTION 

The Applied Mathematics Division at the David Taylor Research Center (DTRC) has begun 
developing methods using finite elements with NASTRAN to solve problems involving 
electromagnetic waves propagating in various media or scattered by objects in the field. This 
paper reports work supported by the Office of Naval Technology Exploratory Development 
Program, DTRC Project Manager, Dr. Bruce Hood. 

The fundamental equations governing the propagation of electromagnetic waves are the 
Maxwell’s equations. For many applications, the electric and magnetic field components 
satisfy the linear, damped wave or Helmholtz equation. While there are six field components 
in electromagnetic problems, for time harmonic fields, only three are independent. The 
equations are, therefore, similar to the Navier’s equations governing an elastic solid. 
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In this paper, the several ways of formulating Maxwell’s equations are presented. 
Formulations involving field vector components are compared with formulations involving 
potential functions. It is shown in the next section that it is possible to form analogies 
between Maxwell’s equations and the Navier’s equations. Standard finite element codes 
which solve the equations of elasticity (such as the NASTRAN code), therefore, with 
appropriate choices of material properties and boundary conditions, can be used to solve 
problems in electromagnetics. 

In this paper, several example problems in electromagnetism are solved using elastic analogies 
and the NASTRAN finite element code. Examples of interest in low frequency applications 
and high frequency radar cross section applications are presented. The examples are all two 
dimensional; however, the analogies and the ability to solve electromagnetism problems with 
NASTRAN are not limited to two dimensional applications. All the applications use the 
IS2D8 element; however, any of the solid elements can be employed. 

The problems of modeling point dipoles in both conducting and nonconducting media are 
studied in this paper. The accurate modeling of dipole sources is critical for applications 
involving electromagnetic waves. The results of these problems are compared with available 
analytic solutions. An important result is determining the mesh requirements needed to 
establish the dipole field accurately. 

The mesh characteristics required for dipole modeling are employed for the study of the fields 
generated from a point dipole source located in sea water (which is a conducting, attenuating 
medium). Two frequencies representing the extremes of low frequency applications are 
presented. Of special interest is the study of the effect of a layer of ice on the solutions. 
While the example presented is somewhat idealized and limited, it should demonstrate to the 
reader the methodology required for the solution of low frequency problems. 

Another example problem is the scattering of a plane wave by a conducting object. The 
problem of a circular cylinder in a plane wave field is solved and compared with the analytic 
solution. Excellent agreement is demonstrated using the IS2D8 element. Convergence of 
this element is quite superior to linear elements documented elsewhere. For this problem, 
both the electric field vector and the magnetic field vector equations are solved. This is 
analogous to the “sound soft” and “sound hard” scattering problems in acoustics. 

Finally, the concluding section of this paper discusses areas where further development is 
required to solve some difficult, three dimensional problems. There is great potential in using 
standard finite element codes for the routine solution of electromagnetics problems. 


MAXWELL’S EQUATIONS - FIELD STRENGTH AND POTENTIAL FORMULATION 

Media in which electromagnetic waves travel often exhibit the properties of linearity, isotropy 
and homogeneity. This type of medium is called a linear, isotropic, and homogeneous (LIH) 
medium. For these media, the electric displacement vector, D, the magnetic field intensity, 
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H, the electric field intensity, E, the magnetic induction vector, B and the free current 
density, J are linearly related by the equations 

D = eE; H = B//i; J = crE (1) 


where e is the permittivity , /i is the permeability and a is the conductivity of the material. If 
we restrict our discussion to time harmonic fields at a single, arbitrary frequency, the field 
vectors E and H take the form 


E = E° exp( ioj t) 
H = H° exp( iwt) 


(2) 


where oj is the frequency and i = \/-T. For time harmonic fields in an LIH medium, the 
governing equations for the spatial field variations can be derived from the Maxwell’s 
equations and are given by [l] 
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where 





( 4 ) 


is the complex permittivity, c is the speed of light in the medium and e 0 is the permittivity of 
free space. The subscript, r, is not summed as is standard in the literature of 
electromagnetism [l] while the subscript i and j are summed in the standard cartesian tensor 
notation These equations are damped wave equations. 

The governing equations given in (3) are the general equations to be solved for any problem 
in electromagnetics involving LIH media. It is important to note that these equations are 
uncoupled. The boundary conditions, however, may involve combinations of the field 
variables. The total problem, therefore, may be strongly coupled. This system, (3), 
represents six partial differential equations in the three components of E and of H. These 
variables, however, are not all independent. For time harmonic applications, only three of 
these equations are independent. For two dimensional time harmonic problems, only two of 
the components are independent. It is important, therefore, to insure that the problem under 
investigation is well posed. In practice, three dimensional electromagnetic problems are 
solved by solving for either E or H and calculating the other from the Maxwell equations. 
For special applications, one could choose two components of one field and one of the other. 
It is important to choose primary unknowns which are consistent with the available boundary 
conditions. 
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An alternative approach to the formulation of electromagnetic problems is to introduce 
potential quantities and derive governing equations for them from the Maxwell equations. A 
vector potential, A and a scalar potential <j> are introduced by the relations 


B = V X A 


E = - y 4> - 


dA 

dt 


( 5 ) 


The potentials are not independent. They can be related to each other using the Lorentz 
gauge condition given by 


v A+£ "^- ”° 


(6) 


This is not the only possible gauge condition relating these quantities; however, it is the most 
widely employed [2]. 

If the relations in (5) are substituted into the Maxwell equations and the vector potential, A, 
is assumed to have a harmonic time variation given by 

A = A°exp(zcot) (7) 


then the governing equation for A is 
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This system, (8), is also a damped wave-like equation with “shear coupling”. The mixed 
partial derivative term comes from the conductance property of the medium. For 
nonconducting media, this equation reduces to an undamped Helmholtz equation. 


NAVIER’S EQUATIONS AND ELASTIC ANALOGIES 

For elastic bodies with time harmonic displacement response, the displacement vector for the 
steady state forced response (at frequency, a >) of the domain is given by 

u = u°exp(zo;t) (10) 
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The governing equation for a damped, isotropic elastic media can be written as [3] 
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where G is the shear modulus, X is the Lame constant, p is the mass density and b is the 
damping coefficient. An alternative form is 
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where 
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(13) 


are complex material parameters. These are the equations which are solved by finite element 
codes designed for the solution of forced, harmonic elastic systems (such as NASTRAN). 

It is desired to draw an analogy between the Navier equations and the Maxwell equations. 
This has been discussed in the literature previously for the scalar Helmholtz equation [4]. 
Following this approach, introduce the relation between Young’s modulus, Y and the shear 
modulus G 


Y = aG = 2(l+i/) G 


(14) 


where the Poisson’s ratio, u, is 


v = 


a 

2~ 


(15) 


If the parameter a is chosen large enough so that 
a+1 ~ a 


(16) 


then 

H, « 0 (17) 

The Navier equations, under this choice of a, reduce to the Maxwell equations of (3). For 
most computers, a value of a = 10 20 is usually sufficient [4]. The shear modulus, G, can be 
chosen arbitrarily. 

If the problem of interest is two dimensional, the Navier equations must be reduced to the 
equations of either plane stress or plane strain. For plane stress, introduce the parameter, j3, 
in the relation 
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Y = fiG 


(18) 


where 


v 



(19) 


If is chosen so that 

/3+1 f* 1 


(20) 


then, for the case of plane stress, [ 4) 
H x « 0 


(21) 


For scalar field problems on most computers, the choice of /? = 10 5 is sufficient [4]. The 
shear modulus, G, can still be chosen arbitrarily. 


For either the two dimensional plane stress analogy or the three dimensional analogy, the 
complex electromagnetic material properties are related to the elastic properties through the 
equation 


He 
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( 22 ) 


The full Maxwell equations for an arbitrary LIH medium (two or three dimensional) can be 
solved by any finite element code which solves the Navier equations if the material properties 
are chosen appropriately and if the boundary conditions can be related to the applied forces 
and displacements. Boundary conditions will be discussed more rigorously in a later section. 

An analogy can be formulated for the magnetic vector potential if the medium is 
nonconducting. In this case, the procedure is identical to the previous discussion as the 
mixed derivative terms do not appear. If the material is conducting, an analogy can be made 
if 

H x = 7 ; H 2 = k (23) 


This is possible if the elastic constants are complex. Since most structural codes do not 
permit complex material constants, the implementation would prove difficult. If one 
examines these material analogies, however, it can be seen that the required complex stiffness 
matrix can be formed by the sum of two real stiffness matrices multiplied by complex 
coefficients. NASTRAN can accomplish this by using DMAP instructions. The imaginary 
part of the stiffness matrix can be calculated in an analogy where 
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ojjjcr = H x ; 


(24) 


li o c 2 


This matrix can be saved, multiplied by - i } and added to the stiffness matrix for a problem 
with 

= 0; H 2 = (25) 

c 

The new stiffness matrix will be the required matrix. 


PIECEWISE HOMOGENEOUS MEDIA 


In many applications, it is necessary to describe the electromagnetic fields which pass from 
one medium to another. Such problems are piecewise homogeneous. For problems with 
only dielectric materials (no conducting materials), this can be done by insuring all elements 
contain only one material and using different element material properties for the different 
media. The procedure is identical to solving problems where the density or elastic modulus 
varies from element to element. 


For conducting media, the “viscous” damping coefficient needs to vary from element to 
element. This is not possible directly with the NASTRAN code. It is possible, through the 
use of DMAP statements, to simulate this with two matrix formulation runs. Form the mass 
matrix for the model with a mass density given by 
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(26) 


The mass density in each element can be different representing the different conducting 
media. It is desired to form a damping matrix, B, which has the element damping coefficients 
given by 

Gum 

b = — produces B (27) 

C «o 

This is accomplished if 

B = M, (28) 


The first matrix formulation produces the M, matrix. This can be written out and read in 
(using OUTPUT2/INPUT2) as the damping matrix in run two which now uses a mass density 


p = ^ - produces M 


(29) 


to form the true analogous mass matrix. Each element can have different permittivities and 
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permeabilities. The only requirement is that the material parameters are constant within an 
element. 

Using the approach summarized above, multiple conducting media can be modeled with 
NASTRAN. An example of a layered media problem is presented in the following section. 
When this procedure is applied in NASTRAN, it is necessary to add a single damper element 
to the model with a zero damping coefficient. This will signal NASTRAN that the problem is 
fully complex and that the complex solver is required. Reading in the damping matrix is not 
sufficient for NASTRAN to choose the complex solver. If other dampers are present in the 
model, this is unnecessary. 


TWO DIMENSIONAL EXAMPLE PROBLEMS 

Several example problems are solved in this section to demonstrate the use of the analogies 
described previously. The problems presented range from very low frequency examples (at 1 
Hz) to high frequency scattering examples (at 3 x 10 8 Hz). All the models employ the 8 
node, quadratic, isoparametric quadrilateral element (IS2D8). These elements perform well 
for a variety of problems and yield accurate results for the problems with available analytic 
solutions. 


EXAMPLE 1: A DIPOLE SOURCE IN FREE SPACE 

As the first problem, the field produced by a two-dimensional point dipole in free space was 
computed to explore the use of analogies with NASTRAN, Information gained by computing 
the fields for this case will also be useful if fields in layered media need to be computed with 
a dipole source located in air. The finite element mesh used is shown in increasing detail in 
Figures 1, 2, and 3. Similar mesh configurations were used with two sizes of elements. For 
the larger elements the overall dimensions of the mesh (Figure 1) are 6T0 8 by 6*10 8 m. 
Thus each of the larger square elements in Figure 1 are 10 8 by 10 8 m, and there are three of 
these elements for each wavelength. The overall dimensions of the smaller mesh are 310 8 by 
310 8 m. For this mesh there are six elements per wavelength. The relative dimensions of all 
elements in the two meshes are equal, so each is portrayed by the figures. The radial mesh in 
the lower left corner of Figure 1, which is graded down to ever smaller elements, contains the 
dipole source. In this section of the mesh, which has dimensions 10 8 by 10 8 m in the larger 
mesh (510 7 by 5*10 7 m in the smaller), the elements are much smaller and the only 
consideration on the element size is to keep the aspect ratios within reasonable bounds (less 
than 1:8). 

Boundary conditions are applied to the model to provide for wave absorption at the outer 
boundary, to apply symmetry conditions on the axes of symmetry, and to model the dipole 
source. The dipole boundary conditions are applied along the small circular boundary in the 
lower left-hand corner of Figure 3. Along the outer boundaries (upper and right sides in 
Figure 1), plane wave absorbing boundary conditions in the form of dashpots were applied. 
Along the axes of symmetry (lower and left sides in Figures 1, 2, and 3), symmetric 
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boundary conditions were applied. The dipole source is modeled by imposing enforced values 
of the electric field for a dipole on the circular boundary sector in the lower left corner of 
Figure 3. The radius of this sector is 0.1 m. The complete solution for a two dimensional 
dipole can be found in [5] 

The electric fields computed for this problem were compared with analytic solutions [5], and 
both meshes were found to produce reasonably accurate values. The amplitude and phase of 
the solution for the larger mesh are shown in Figures 4 and 5. and the amplitude and phase 
of the near-field solution are shown in Figures 6 and 7. The values plotted are the z- 
component of the electric field along a radial line 45 degrees from the lower axis. For the 
larger mesh, the error in the large square elements was on the order of 5 percent, and in the 
radial block the error was of the order of 1 percent. For this model, the region containing the 
radial elements is considered to be the region of interest, and the outside region is included 
only to model several wavelengths to provide for suitable wave absorbing boundaries. 
Therefore, in the region of interest, very good results were obtained. 

For the larger mesh, two wave lengths were modeled before the absorbing boundary 
conditions are applied, and for the smaller mesh only one wave length was modeled. 
Decreasing the number of wave lengths modeled inside the boundary increased the error in 
the radial elements, the region of interest. The change in mesh size resulted in errors of 4 
percent in both the square and radial elements in the smaller mesh. At the same time, the 
increase in the number of elements per wave length in the outer region slightly increased the 
accuracy there. 


EXAMPLE 2: A DIPOLE SOURCE IN SEA WATER 

Computing the field due to a dipole source in sea water was the first application to modeling 
electromagnetic fields in a conducting medium. As with the preceding problem, the region 
containing seawater was assumed to have infinite extent, so that comparisons could be made 
to an analytic solution [5], Since for frequencies near one Hertz, sea water is a good 
conductor, the electromagnetic wave length, equal to 1581 m, is considerably shorter than 
3T0 8 in free space. Therefore, the region modeled for this problem was correspondingly 
smaller than the region for the preceding problem. The finite element mesh used is shown in 
increasing detail in Figures 8, 9, and 10. The outer dimensions of this mesh are 5000 by 
5000 m. In the outer region of Figure 8, the larger square elements are 250 m on aside, and 
the smaller square elements on the left side are 125 m on a side. The elements on the left 
were made smaller because this same mesh was to be used as part of the layered media 
problem, and the use of various element sizes allowed checking the performance of 
transitions from smaller to larger elements. The radial mesh in the lower left corner of 
Figure 8 contains the dipole source which is too small to be seen in this figure, but can be 
seen in Figure 10. This section of the mesh has dimensions 125 by 125 m. Here again the 
elements are much smaller, and the principal consideration is to keep the aspect ratios 
reasonable. 
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Again, boundary conditions are applied to provide for wave absorption on the outer 
boundary, to apply symmetry conditions on the axes of symmetry, and to model the dipole 
source. The dipole boundary conditions are applied along the small circular boundary in the 
lower leftrhand corner of Figure 10. Along the outer boundaries (upper and right sides in 
Figure 8), plane wave absorbing dashpots were applied. Along the axes of symmetry (lower 
and left sides in Figures 8, 9, and 10), symmetric boundary conditions were applied. The 
dipole source is modeled by imposing enforced values of the electric field for a dipole on the 
circular boundary sector in the lower left corner of Figure 10. The radius of this sector is also 
0.1 m. 

The electric fields computed for this problem were compared with analytic solutions, and were 
found to produce accurate values. The amplitude and phase of the solution along the 
horizontal axis of symmetry are shown in Figures 11 and 12. The solution phase is plotted 
between -180 degrees and 180 degrees, therefore, an apparent discontinuity arises at radii at 
which the phase decreases past -180 degrees. The error in the solution was on the order of 1 
percent everywhere. Again for this model, the region containing the radial elements is 
considered to be the region of interest. The outside region is included only to model enough 
of the medium to provide for wave absorbing boundaries. Therefore very good results were 
obtained in the region of interest. The amplitude and phase of the solution along the vertical 
axis of symmetry are shown in Figures 13 and 14. Excellent agreement with the analytic 
solution is demonstrated in this direction also. 


EXAMPLE 3: MODELING A DIPOLE SOURCE IN A FINITE DEPTH OF SEA WATER 

The problem under consideration is a dipole source located in a finite depth of sea water. 
The current modeling is limited to a two-dimensional line dipole. The general problem under 
consideration is shown in Figure 15. A two-dimensional dipole is located at a distance A 
beneath the surface of the sea water. The total depth of the sea water is H. The sea water 
may be covered by a layer of ice of thickness D. The air on top and the mud beneath the sea 
water are assumed to be infinite. The problem currently modeled assumes a sea depth of 250 
m. The dipole source is located 125 m beneath the surface of the sea water. Models have 
been developed for sea water without ice and for sea water covered by 1 m of ice. The total 
mesh for the problem under consideration for a dipole source radiating at 1 hertz is shown in 
Figure 16. 

The sea water is modeled for a total of 7500 m (approximately 30 skin depths) and then is 
terminated by a plane wave radiation boundary condition. The mud is modeled out to 16000 
by 16000 m (which is approximately 20 skin depths). In the model shown, the air is also 
modeled out to 16000 by 16000 m. Both media are terminated by plane wave radiation 
boundary conditions. Since mud and sea water are attenuating media, the radiation boundary 
condition assumption is not expected to significantly influence the solution (this has been 
demonstrated for the case of a line dipole in an infinite region of sea water as discussed 
previously). For the air, however, it is often required to model a region on the order of 
several wave lengths. For air (at 1 hertz excitation), this corresponds to approximately 
300,000,000 m. Results in air, however, are only of interest for distances less than 10,000 m 
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from the source. The air, therefore, was modeled as far as the mud (for geometric 
symmetry) and absorbing boundary conditions were applied at the edge of the mesh. 

Figure 17 shows a blowup of the entire sea water region. The transitioning mesh in the air 
and the sea water is shown. For the dipole in sea water, a mesh dimension of 125 by 125 m 
was demonstrated to predict accurate results. This is the dimension of the elements in the 
sea water as shown. In the air and mud, the elements are allowed to expand in a consistent 
manner to a final dimension of 1000 by 1000 m. The larger elements are permissible since 
the wavelength and attenuation distance in mud are larger than in seawater (the wavelength 
in mud is in the order of 5000 m and the skin depth in mud is on the order of 796 m). The 
transitioning is developed to insure that the element aspect ratios and interior angles remain 
within acceptable limits. Near the dipole source, a radially expanding mesh is employed as in 
the previous example. This mesh is sufficient to establish the near source field accurately. 
As in the previous example, the dipole is modeled as a small circular ring of nodes. On that 
ring of nodes, the analytic solution for a line dipole in an infinite medium of sea water is 
applied as a boundary condition. The model assumes, therefore, that close enough to the 
dipole, the ice, air and mud will have a negligible effect on the field variable solutions. The 
required inner mesh dimension will be determined by a convergence study. This parameter 
will be dependent upon the location of the source relative to the boundaries and the 

frequency of the source. The model described was modified to allow for a 1 m layer of ice. 

The resulting mesh is the same as the previous one except that between the sea water - air 
interface is a layer of elements 1 m thick which represent the ice. Since this dimension is 
small relative to the domain modeled, it is observable only on a blowup of the mesh. 

Figures 18 and 19 show the amplitude and phase of the electric field component, E z , along 
the midline of the sea water. Solutions with and without ice are shown. A decaying field is 
observed with a characteristic knee in the solution. This occurs near the point where the 
phase crosses the zero line. This phenomena has also been observed experimentally [6j. The 
dropoff in the phase near the tail of the plot is probably due to the dash pot boundary 
condition. Figures 20 and 21 show the amplitude and phase of the E z component along the 
surface of the sea water. Qualitatively, the solution is similar to the midline solution. The 
amplitude does not, however, drop off as rapidly and the phase is shifted to a larger mean. It 
is interesting to note that at this low frequency the ice has negligible effect on the solution. 

The same problem was studied for a higher frequency source at 1000 hz. The mesh 
employed is shown in Figure 22. The seawater region is modeled for 250 m by 250 m. This 
corresponds to about 50 skin depths. The mud is modeled for an additional 500 m 
corresponding to 20 skin depths in mud. The air is modeled out to 2000 additional meters. 
On all exterior boundaries, the dashpot absorbing conditions are employed. 

The amplititude and phase of the E z component are shown in Figures 23 and 24 along the 

midline of the sea water. Again, solutions with and without ice are shown. The solution has 

a typical decaying amplitude with a sawtoothed phase characteristic. This is similar to the 
solution for a dipole in a conducting medium. It is interesting to note, however, that while 
the wavelength corresponds to the wavelength of the media, the decay is slower than for a 
dipole in infinite sea water. The amplitude is receiving significant contribution from the 
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surface waves along the sea surface. Along the sea midline, little influence of the ice can be 
seen. 

Figures 25 and 26 show the amplitude and phase of the E z component along the surface of 
the sea water (with and without ice). The ice clearly has a significant influence on this 
solution. The amplitude without ice follows the amplitude with ice for about 10 skin depths 
of the sea water. The two solutions then change and the amplitude with ice is larger. Even 
though ice has a small conductivity (10 -5 mhos), it acts as a wave guide keeping the surface 
wave of larger amplitude than without the ice. The phase, however, shows little difference 
with and without ice. The solutions are qualitatively similar with the exception that the ice 
guides the surface wave. Note that the deviation of the phase and the slight increase in 
amplitude toward the end of the plots is probably due to reflections from the dashpots. 


EXAMPLE 4: PLANE WAVE SCATTERING FROM A RIGHT CIRCULAR CYLINDER 

As a final example, consider the scattering of an incident plane wave by an infinite, perfectly 
conducting circular cylinder. The boundary condition on the cylinder is that the longitudinal 
component of the electric field must vanish on the surface of the cylinder and that the 
longitudinal component of the magnetic field must be normal to the surface. If the governing 
equations for the scattered wave only are considered, the boundary conditions for the 
scattered wave must remove the E z component of the plane wave at the cylinder surface. In 
addition, the normal derivative of the H z component of the plane wave must vanish at the 
surface of the cylinder. At infinity, the scattered wave must vanish. The problem considered 
is for aim cylinder with an incident plane wave of 1 m wave length (the frequency, 
therefore, is 3 x 10 8 Hz). 

The first mesh attempted employed eight elements in the azimuthal direction and quarter 
wavelength dimension in the radial direction. As is shown subsequently, this mesh 
performed adequately for the longitudinal component of the magnetic field but was not 
sufficient to accurately solve the longitudinal electric field problem. The mesh employed for 
the longitudinal electric field component is shown in Figure 27. The mesh was generated 
using the IDEAS [7] package. The design criterion was to generate a mesh as close to 
uniform in dimension as possible with an element size equal to one quarter of the incident 
wave length. The performance of this mesh was superior to that of a mesh with fixed radial 
dimensions of one quarter of a wavelength and aspect ratios within 1 to 5. The zero field 
condition was modeled with absorbing dashpots. For this case, a cylindrical wave condition 
would be superior due to the geometry of the problem. This was compared with the simple 
dashpot condition for the magnetic field solution. 

Figure 28 is a plot of the normalized amplitude of the scattered electric field intensity on the 
forward scattering side of the cylinder. The overall agreement is quite good. The maximum 
error is less than 5% compared with the analytic solution. Very near the cylinder, however, 
the largest deviation is observed. Indeed it is only in this region where the error is larger 
than 1%. Figure 29 is a blowup of this region. The analytic solution flattens near the cylinder 
while the finite element solution demonstrates a sharp dip. There is significant ripple in the 
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solution in this region which may indicate reflection problems between the cylinder and the 
dashpot. A blowup of the analytic solution is shown in Figure 30. The region which 
appeared flat in Figure 29 has a slight dip as demonstrated in Figure 30. The finite element 
solution exaggerates this dip. Since the elements are on the order of 0.25 m (one quarter of a 
wavelength), it is evident from Figure 30 that this mesh density would be insufficient to 
totally reproduce this phenomena. Overall, however, the solution is quite good. 

Figure 31 shows a plot of the phase of the forward scattered field. The finite element results 
are almost identical to the analytic solution. This demonstrates that while small amplitude 
errors may be introduced into the solution, the general character of the waves are accurately 
predicted by the finite element solution. Figure 32 shows the normalized amplitude of the 
scattered electric field on the back side of the cylinder. The finite element results agree very 
well with the analytic solution. All errors are bounded by 1 % even near the dashpot 
boundary condition. Figure 33 shows a plot of the phase of the electric field on the 
backscattering side of the cylinder. Again, excellent agreement is seen. The sawtoothed 
phase characteristic is accurately predicted and the ramping behavior is accurate. 

For this example, the longitudinal component of the magnetic field vector was also resolved. 
The longitudinal components of the E and H fields are the only independent components for 
the two dimensional applications. Figure 33 shows the amplitude of the H z component as a 
function of distance away from the cylinder along the back scattering side. The finite element 
solution is only negligibly different from the analytic solution. Figure 34 shows the amplitude 
along a radial line at 112.5 degrees from the incident wave. This represents the worst case 
and yet the two solutions agree quite well. It should be noted that the total amplitude along 
this line is quite small. It is remarkable that the solution is this accurate. In addition, the 
finite element results quite accurately capture the spiked dip in the solution even though only 
four elements per wavelength were employed. 

An important observation is that accurate solutions were generated with approximately four 
IS2D8 elements per wavelength. This problem has been solved previously with linear 
quadrilateral elements [8]. In that study, ten elements per wavelength were required 
necessitating a significantly greater number of degrees of freedom to achieve an accurate 
solution. 


CONCLUDING REMARKS 

Maxwell’s equations were solved for a variety of example problems in two dimensions. An 
interesting outcome of the low frequency examples was the ability to predict the wave guide 
effects of the ice in the layered media problem and the knee response in the amplitude. In 
addition, this problem demonstrated that relatively complicated problems can be solved by 
routine methodology. This conclusion should also hold for three dimensional applications. 

The scattering example demonstrates the ability to handle high frequency applications. The 
conclusion that only four IS2D8 elements are required per wavelength indicates that 
considerable economy should be realized by using quadratic elements for harmonic response 
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applications. This conclusion should be valid for structural and acoustic applications as well as 
for electromagnetic applications. 

Since the two dimensional problems exhibit totally uncoupled boundary conditions, the 
solution of two dimensional problems in electromagnetics is the same as solution of the scalar 
wave equation. In three dimensions, this is not the case. The boundary conditions 
encountered are often coupled. This poses a problem for certain situations. 

A common three dimensional boundary condition is the perfect conductor condition of zero 
tangential E and normal H. For high frequency applications, this is not a problem because 
several skin depths of the conductor can be modeled easily since this dimension will be small 
relative to the conductor’s size. The conductor will damp out and absorb the waves 
appropriately. This indicates that radar cross section problems in three dimensions can be 
handled by elastic analogies. 

For low frequency applications, the presence of a conductor is not as easy to deal with since 
the skin depth is often large relative to the size of the conductor. The vanishing of the 
tangential E field can be handled by multipoint constraints (MPCs). The vanishing of the 
normal H field is not as trivially solved. Methodologies for enforcing this condition are under 
investigation. It may be possible to extend the concept of MPCs to include linear 
combinations of first partial derivatives. This would solve the problem. 

The other major problem to be addressed is the fact that many electromagnetic problems are 
exterior problems. They involve either extremely large or infinite domains. The solution is 
of interest, however, only in a small domain. It is necessary, therefore, to reduce the 
modeled domain and to implement a boundary condition which accounts for the remaining 
media. In this paper, the simple plane wave condition was employed. While this works, 
often large domains must be modeled. Other conditions have been explored; however, 
additional research is required. Infinite elements (employed for some limited scalar 
applications [9]) hold promise. These are currently being investigated also. 

The remainder of the boundary conditions encountered in most applications can be handled 
trivially with elastic finite element codes like NASTRAN. This paper has demonstrated the 
ability to handle two dimensional problems and has provided the formulation for three 
dimensional problems. When absorbing boundary conditions become available and the zero 
normality condition for the H field is developed, it will be possible to solve virtually all 
problems in electromagnetics (which adhere to the assumptions of Maxwell’s equations) with 
elastic finite element codes like NASTRAN. 
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Fig. 5. PHASE FOR POINT DIPOLE IN FREE SPACE 



Fig. 6. NEAR FIELD AMPLITUDE FOR POINT DIPOLE IN FREE SPACE 
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Fig. 7. NEAR FIELD PHASE FOR POINT DIPOLE IN FREE SPACE 



Fig. 8. FINITE ELEMENT MODEL FOR POINT DIPOLE IN SEA WATER 





Fig. 9. RADIAL DIPOLE MODEL IN SEA WATER 



233 






AMPLITUDE Of E- 



235 



AIR 




Fig. 16 . FINITE ELEMENT MODEL FOR DIPOLE IN SEA WATER WITH AIR, MUD 
AND ICE 
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Fig. 22. FINITE ELEMENT MESH FOR 1000 HZ. SOURCE 
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Fig. 23. AMPLITUDE OF E-Z ALONG SEA MIDLINE - 1000 HZ. 



Fig. 24. PHASE OF E-Z ALONG SEA MIDLINE - 1000 HZ. 
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Fig. 26. PHASE OF E-Z ALONG SEA SURFACE - 1000 HZ. 
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Fig. 27. FINITE ELEMENT MESH FOR SCATTERING ABOUT CIRCULAR CYLINDER 



Fig. 28. AMPLITUDE OF THE FORWARD SCATTERED E-Z WAVE 
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Fig. 29. BLOWUP OF THE SCATTER IN THE E-Z AMPLITUDE SOLUTION 



Fig. 30. NEAR CYLINDER BEHAVIOR IN THE ANALYTIC SOLUTION 
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Fig. 31. PHASE OF THE FORWARD SCATTERED E-Z WAVE 



Fig. 32. AMPLITUDE OF THE BACICSCATTERED E-Z WAVE 
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Fig. 35. WORST CASE FOR THE AMPLITUDE OF THE H-Z WAVE 
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